In [1]:
#import relevant libraries 
import numpy as np
import pandas as pd
from scipy.optimize import root_scalar
from scipy import interpolate
import plotly.graph_objects as go
import plotly.figure_factory as ff
In [2]:
# create functions for duct curvature and an array of x values to use in graphing 
x = np.linspace(0, 1, 500)
alpha = (x/2 - 0.3)**2 + 0.05 
beta = -alpha

# define given variables and constants 
p_0i = 140 # kPa
T_0i = 290 # K
R = 287 # J/kgK
c_p = 1006 # J/kgK
g = 1.4
shock_x = 0.8 
throat_x = 0.6

# assume that flow is choked so the throat Mach number can be defined as 1
M_t = 1

Information to keep in mind about the flow: The flow is given to be in a frictionless duct and is a perfect gas The flow before and after the shock will be isentropic and therefore isentropic flow relations can be used \ Across a normal shock the stagnation pressure will drop and stagnation temperature will remain the same \ Due to conservation of mass the mass flow rate is constant through out the duct

The curvature of the top of the nozzle is given by $ \alpha(x) $ and the curvature of the bottom is the negative of the top: $$ \alpha(x) = \left(\frac{x}{2}-0.3\right)^{2} + 0.05 = -\beta(x) $$

$x$ is given to be $ x \in [0, 1] $

In [3]:
# Define a function that solves for an area per unit width (can also solve for different areas such as circular)
def area(x):
    area = 2 * ((x/2 - 0.3)**2 + 0.05)
    return area

To find the Mach number in the nozzle before the shock we can use the isentropic area ratio mach number relation: $$ \require{color}\frac{A}{A^{*}} = \frac{\left(1 + \frac{\gamma+1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^{2}\right)^{\frac{\gamma+1}{2(\gamma-1)}}}{{\color[rgb]{0.041893,0.355669,0.727621}M}} \left(\frac{\gamma +1}{2}\right)^{-\frac{\gamma +1}{2(\gamma -1)}} $$

There is no analytical solution to find the Mach number from the area ratio so another method must be used, in this notebook optimization is used

In [4]:
# Define the area at the throat, which is needed in order to solve for Mach number
A_t = area(throat_x)
In [5]:
# Define a mach function to calculate the roots of the area ratio mach relation set to zero
def mach_function(M, x_pos):
    A = area(x_pos)
    return np.sqrt((1/M**2)*((2/(g+1))*(1+((g-1)/2)*M**2))**((g+1)/(g-1))) - A/A_t
In [6]:
# Define the Mach number at the shock position, min Mach value at the entrance, and Mach number after shock
M_s_root = root_scalar(mach_function, args = (shock_x), x0 = 1, x1 = 2.5)
M_s = M_s_root.root

M_min_root = root_scalar(mach_function, args = (0), x0 = 0.2, x1 = 1)
M_min = M_min_root.root

M_as = np.sqrt(((g-1)*M_s**2 + 2)/(2*g*M_s**(2)-(g-1)))

To find the Mach number after the shock we need the mass flow rate: $$ \require{color}\dot{{\color[rgb]{0.501961,0.250953,0.010028}m}} = \frac{A{\color[rgb]{0.315209,0.728565,0.037706}p}_{0}}{\sqrt{{\color[rgb]{0.121820,0.954406,0.966585}T}_{0}}} \sqrt{\frac{\gamma}{{\color[rgb]{0.986252,0.007236,0.027423}R}}} {\color[rgb]{0.041893,0.355669,0.727621}M} \left(1+\frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^{2}\right)^{-\frac{\gamma+1}{2(\gamma -1)}} $$

This equation also cannot be solved analytically so the same process as before can be used \ Since the mass flow rate must remain constant to obey conservation of mass and the stagnation temperature also remains constant, the area and the stagnation pressure are the only changing variables

The stagnation pressure will remain constant at the inlet value until the flow reaches the shock Because the stagnation pressure after the shock only depends on the Mach number at the shock and the inlet stagnation pressure we can go ahead and define a function to calculate stagnation pressure based on the normal shock relation: $$ \require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}_{0s} = {\color[rgb]{0.315209,0.728565,0.037706}p}_{0i} \left(\frac{(\gamma+1){\color[rgb]{0.041893,0.355669,0.727621}M}_{s}^{2}}{(\gamma-1){\color[rgb]{0.041893,0.355669,0.727621}M}_{s}^{2}-2}\right)^{\frac{\gamma}{\gamma-1}} \left(\frac{\gamma+1}{2 \gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_{s}^{2} - (\gamma -1)}\right)^{\frac{1}{\gamma-1}} $$ Where $\require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_{s} $ is the value of the Mach number at the shock

In [7]:
# Define a function to calculate the stagnation pressure before and after the shock -- based on shock relations
def stag_pressure(x_pos):
    if x_pos < shock_x:
        stag_pressure = p_0i       
    else:
        stag_pressure = p_0i * ((((g + 1) * M_s**2) / ((g - 1) * M_s**2 +2))**(g/(g - 1))) * ((g + 1) / (2*g* M_s**2 - (g -1)))**(1/(g-1))
    return stag_pressure
In [8]:
# Define the mass flow rate at the throat where M=1 for ease, this number will remain constant
m_dot = ((A_t)*(p_0i*1000)/(np.sqrt(T_0i)))*np.sqrt(g/R)*((g+1)/2)**(-(g+1)/(2*(g-1))) 
#p0 is multipled by 1000 because it is given in kPa and we need units of Pa

# Define a function that finds the roots of the mass flow rate equation set to zero
def mach_function_as(M, x_pos):
    A = area(x_pos)
    p_0 = stag_pressure(x_pos)*1000
    return (A*p_0/np.sqrt(T_0i))*(np.sqrt(g/R))*M*(1+((g-1)/2)*M**2)**(-(g+1)/(2*(g-1))) - m_dot
In [9]:
# Define a Function to calculate the Mach value at an x value, x0 and x1 are lower and upper limits on the Mach number that
# are based on location in the nozzle
def mach_number(x_pos):
    if x_pos <= 0.5:
        M_root = root_scalar(mach_function, args = (x_pos), x0 = M_min, x1 = 0.8)
        M = M_root.root
    elif 0.5 < x_pos <= throat_x: 
        M_root = root_scalar(mach_function, args = (x_pos), x0 = M_min, x1 = 1)
        M = M_root.root
    elif x_pos < shock_x:
        M_root = root_scalar(mach_function, args = (x_pos), x0 = 1, x1 = M_s)
        M = M_root.root
    else: 
        M_root = root_scalar(mach_function_as, args = (x_pos), x0 = M_as, x1 = 0.2)
        M = M_root.root
    return M

After you have a function that defines the stagnation pressure and Mach number everywhere the static pressure can be calculated using isentropic flow relations, as flow is only momentarily nonisentropic right at the shock so there is a discontinuity: $$ \require{color}{\color[rgb]{0.315209,0.728565,0.037706}p} = {\color[rgb]{0.315209,0.728565,0.037706}p}_{0} \left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^{2}\right)^{-\frac{\gamma}{\gamma-1}} $$

In [10]:
# Define a function to calculate the static pressure -- based on isentropic flow relations 
def static_pressure(x_pos):
    M = mach_number(x_pos)
    p_0 = stag_pressure(x_pos)
    static_pressure = p_0 * (1 + ((g - 1)/2)*M**2)**(-g/(g-1))
    return static_pressure

To calculate the static temperature isentropic flow relations can also be used, keeping in mind that the stagnation temperature remains constant over a shock: $$ \require{color}{\color[rgb]{0.121820,0.954406,0.966585}T} = {\color[rgb]{0.121820,0.954406,0.966585}T}_{0} \left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^{2}\right)^{-1} $$

Where $ T_{0} $ is the value of the inlet stagnation temperature

In [11]:
# Define a function to calculate the static temperature -- based on isentropic flow relations
def static_temperature(x_pos):
    M = mach_number(x_pos)
    static_temperature = T_0i * (1 + ((g - 1)/2)*M**2)**(-1)
    return static_temperature

Now that the Mach number and static temperature can be calculated at every point, the speed of sound and thus the flow speed at each point can be calculated. To calculate the speed of sound we can use: \begin{equation} \require{color}{\color[rgb]{0.989013,0.435749,0.811750}a} = \sqrt{\gamma {\color[rgb]{0.986252,0.007236,0.027423}R} {\color[rgb]{0.121820,0.954406,0.966585}T}} \end{equation}

To calculate the speed from this equation we can use the relation between Mach number and speed of sound, and solve for speed: \begin{equation} \require{color}{\color[rgb]{0.041893,0.355669,0.727621}M} = \frac{{\color[rgb]{0.059472,0.501943,0.998465}v}}{{\color[rgb]{0.989013,0.435749,0.811750}a}} \end{equation}

In [12]:
# Define a function to calculate the speed of the particle at each position
def speed(x_pos):
    M = mach_number(x_pos)
    T = static_temperature(x_pos)
    a = np.sqrt(g*R*T)
    v = M * a
    return v
In [13]:
# now using the five defined functions from above the values of the five fields of interest can be found at x values
M_values = []
p0_values = []
p_values = []
T_values = []
v_values = []
for x_val in x:
    M_values.append(mach_number(x_val))
    p0_values.append(stag_pressure(x_val))
    p_values.append(static_pressure(x_val))
    T_values.append(static_temperature(x_val))
    v_values.append(speed(x_val))
In [14]:
# using the dataframe we can also create a plot with a dropdown menu between the fields
fig = go.Figure()

fig.add_scatter(x=x, y=M_values, line_color = 'plum')
fig.add_scatter(x=x, y=p0_values, line_color = 'mediumaquamarine', visible=False)
fig.add_scatter(x=x, y=p_values, line_color = 'cornflowerblue', visible=False)
fig.add_scatter(x=x, y=T_values, line_color = 'tomato', visible=False)
fig.add_scatter(x=x, y=v_values, line_color = 'darkmagenta', visible=False)

fig.update_layout(updatemenus = [dict(active=0, buttons = [dict(label='Mach Number', method='update', args=[{'visible':
                  [True, False, False, False, False]}, {'yaxis':{'title':'Mach Number'}}]), dict(label='Stagnation Pressure',
                  method='update', args=[{'visible': [False, True, False, False, False]}, {'yaxis':{'title':
                  'Stagnation Pressure (kPa)'}}]), dict(label='Static Pressure', method='update', args=[{'visible': [False, 
                  False, True, False, False]}, {'yaxis':{'title':'Static Pressure (kPa)'}}]), dict(label='Static Temperature'
                  , method='update', args=[{'visible': [False, False, False, True, False]}, {'yaxis':{'title':
                  'Static Temperature (K)'}}]), dict(label='Velocity', method='update', args=[{'visible':[False, False, 
                  False, False, True]}, {'yaxis':{'title':'Velocity (m/s)'}}])])])

fig.update_yaxes(title_text='Mach Number')
fig.update_xaxes(title_text="x position")
fig.show()

Next we can make contour plots of the fields within the nozzle \ because there is no y variation these plots will be uniform in the y-direction, however when we get to topics such as boundary Layers this will change and will no longer be uniform \ In reality, isentropic flow is an idealized process as is a frictionless duct \ Even in flows with a very low viscosity the viscous forces will never be negligible near a wall and will cause a boundary layer to form thus the field contours would differ along y and leads to the potential for flow separation \ Friction also causes a loss of heat meaning the flow would no longer be adiabatic and causes a drop in enthalpy, therefore reducing efficency as well as the exit Mach number

In [15]:
# create contour maps of the nozzle
fig = go.Figure()

fig.add_contour(z=[M_values, M_values], x=x, contours_coloring = 'heatmap', line_width = 0, colorscale = 'deep',
                reversescale = True, colorbar=dict(title='Mach Number', titleside='right'))
fig.add_contour(z=[p0_values, p0_values], x=x, contours_coloring ='heatmap', line_width = 0, colorscale = 'Teal', 
                visible = False, colorbar=dict(title='Stagnation Pressure (kPa)', titleside='right'))
fig.add_contour(z=[p_values, p_values], x=x, contours_coloring = 'heatmap', line_width = 0, colorscale = 'viridis', 
                visible = False, colorbar=dict(title='Static Pressure (kPa)', titleside='right'))
fig.add_contour(z=[T_values, T_values], x=x, contours_coloring = 'heatmap', line_width = 0, colorscale = 'blackbody', 
                reversescale = True, visible = False, colorbar=dict(title='Static Temperature (K)', titleside='right'))

fig.add_scatter(x=x, y=alpha+0.15, line_color = 'white')
fig.add_scatter(x=x, y=x-x+0.3, fill='tonexty', fillcolor='white', line_color='white')
fig.add_scatter(x=x, y=beta+0.15, fill='tozeroy', fillcolor='white', line_color='white')

fig.update_layout(showlegend=False, updatemenus=[dict(active=0, buttons = [dict(label='Mach Number', method='update', args=
                  [{'visible':[True, False, False, False, True, True, True]}]), dict(label='Stagnation Pressure', 
                  method='update', args=[{'visible':[False, True, False, False, True, True,True]}]), dict(label=
                  'Static Pressure', method='update', args=[{'visible':[False, False, True, False, True, True, True]}]),
                  dict(label='Static Temperature', method='update', args=[{'visible':[False, False, False, True, True, True,
                  True]}])], y=1.2, x=1.1)])

fig.update_yaxes(range=[0, 0.3], showticklabels = False)
fig.update_xaxes(title='x position')

fig.show()

From the velocity function we can visualize how fluid would move through the nozzle. By utilizing the equation $\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v} = \frac{dx}{dt}$ we can find the approximate amount of time between each equally spaced x value within the defined array. \ From time we can find x(t) and linerally interpolate for equally spaced values of t to find approximately how far a particle would move within that time interval

In [16]:
# find the interval of time between each equally spaced x position 
dx = np.diff(x)
vav = []
for i in range(0, 499):
    av = sum(v_values[i: i+1]) / 2
    vav.append(av)
dt = dx / vav
In [17]:
# find the time at each equally spaced x position 
time = []
for j in range(0, 500):
    time.append(sum(dt[0:j]))
In [18]:
# interpolate using the x and time values then use the interpolation function to find x positions at equal intervals of time
f = interpolate.interp1d(x=time, y=x)
tnew = np.linspace(0, max(time), 500)
xnew = f(tnew)
In [19]:
# using the interpolated data we can animate how a particle would move through the nozzle
fig = go.Figure()

fig.add_scatter(x=x, y=np.zeros(500))
fig.add_scatter(x=x, y=np.zeros(500), line_dash = 'dot', line_color = 'lightsteelblue')

fig.add_scatter(x=x, y=alpha, line_color = 'black')
fig.add_scatter(x=x, y=beta, line_color = 'black')
fig.add_vline(x = shock_x, line_dash = 'dash', line_color = 'red', annotation_text = 'shock')
fig.add_vline(x = throat_x, line_dash = 'dash', line_color = 'blue', annotation_text = 'throat')

fig.update_layout(xaxis=dict(range=[0,1], autorange=False), yaxis=dict(range=[-0.15, 0.15], autorange=False), updatemenus=
                  [dict(type='buttons', buttons=[dict(label='play', method='animate', args=[None, {'frame':{'duration':10}}])])])
fig.update(frames = [go.Frame(data=[go.Scatter(x=[xnew[k]], y=[np.zeros(500)[k]], mode='markers', marker=dict(color='teal',
                    size=10))]) for k in range(500)], layout_showlegend=False)

fig.show()

https://www.youtube.com/watch?v=I8ntKdaKxV4 \ This link provides good graphics similar but more involved than mine about how flow moves through a converging diverging nozzle, moving the shock location

In [20]:
# define x and v values for 18 different positions to create arrows
x_18 = np.linspace(0,1,18)
v_18 = []
for x_val in x_18:
    v_18.append(speed(x_val) / 1000)
In [21]:
alpha_2 = (x_18/3 -0.2)**2 + 0.025
beta_2 = -alpha_2
slope_2 = 2*(x_18/3 - 0.2)/3 * abs(np.append(np.sin(np.arctan(np.diff(alpha_2) / np.diff(x_18))), 0))
In [22]:
# create velocity vector field 
fig = ff.create_quiver(x=np.linspace(0,1,18), y=np.zeros(18), u=v_18, v=np.zeros(18), line_color = 'blue')
fig2 = ff.create_quiver(x=x_18, y=alpha_2, u=v_18, v=slope_2, line_color = 'blue')
fig3 = ff.create_quiver(x=x_18, y=beta_2, u=v_18, v=-slope_2, line_color = 'blue')

fig.add_scatter(x=x, y=alpha, line_color = 'black')
fig.add_scatter(x=x, y=beta, line_color='black')
fig.update_layout(xaxis=dict(range=[0,1]), yaxis=dict(range=[-0.15, 0.15]), showlegend = False)

fig.add_traces(data = fig2.data)
fig.add_traces(data = fig3.data)
fig.show()
In [ ]:

all isentropic flow and normal shock equations are taken from NASA